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This paper provides a detailed study of turbulent statistics and scale-by-scale bud- 
gets in turbulent Rayl eigh-Bena r d con vection. It aims at testing the applicability of 
iKolmogorovl l)l94l|) and lBolgianol (|l959^l theories in the case of turbulent convection and 
at improving the understanding of the underlying inertial range scalings, for which a 
general agreement is still lacking. Particular emphasis is laid on anisotropic and inho- 
mogencous effects, which are often observed in turbulent convection between two diffcr- 
entially heated plates. For this p urpose, the SO (3) decomposition of structure functions 
llArad. L'vov fe Procacrialll9996l) and a method of description of inhomogencitics pro- 
posed b\ TDMiaiia^r^ir(f200lJl arc used to derive inhomogeneous and anisotropic gener- 
alizations of Kolmogorov and Yaglom equations applying to Raylcigh-Bcnard convection, 
which can be extended easily to other types of anisotropic and/or inhomogeneous flows. 
The various contributions to these equations are computed in and off the central plane 
of a convection cell using data produced by a Direct Numerical Simulation of turbulent 
Boussinesq convection at Ra = 10 6 and Pr = 1 with aspect ratio A = 5. The analy- 
sis of the isotropic part of the Kolmogorov equation demonstrates that the shape of the 
third-order velocity structure function is significantly influenced by buoyancy forcing and 
large-scale inhomogencitics, while the isotropic part of the mixed third-order structure 
function ((A6>) 2 Ait) appearing in Yaglom equation exhibits a clear scaling exponent 1 in a 
small range of scales. The magnitudes of the various low £ degree anisotropic components 
of the equations are also estimated and are shown to be comparable to their isotropic 
counterparts at moderate to large scales. The analysis of anisotropics notably reveals that 
computing reduced structure functions (structure functions computed at fixed depth for 
correlation vectors r lying in specific planes only) in order to reveal scaling exponents 
predicted by isotropic theories is misleading in the case of fully three-dimensional turbu- 
lence in the bulk of a convection cell, since such quantities involve linear combinations 
of different £ components which are not negligible in the flow. This observation also in- 
dicates that using single points measurements together with the Taylor hypothesis in 
the particular direction of a mean flow to test the predictions of asymptotic dimensional 
isotropic theories of turbulence or to calculate intcrmittency corrections to these theories 
may lead to significant biases for mildly anisotropic three-dimensional flows. A qualitative 
analysis is finally used to show that the influence of buoyancy forcing at scales smaller 
than the Bolgiano scale is likely to remain important up to Ra = 10 9 , thus preventing 
Kolmogorov scalings from showing up in convective flows at lower Rayleigh numbers. 
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1. Introduction 

The quest for inert ial range scal i ng laws in turbulen t conv e ction has been very activ e 
in recent years (e. q . IChilla et all (Il993l) : iBenzi et all l)l994|) ; ICalzavarini et aZ.I (|2002^ : 
IVerzicco fc Camussil (1200.^ : IChing et all 1)2004) ). Their determination is expected to give 
some important insight into the thermal and mechanical processes at work in the flow. 
To this end, turbulent thermal convection is investigated in convection cells heated from 
below using laboratory and numerical experiments, within the framework of the Boussi- 
nesq approximation. Such a flow exhibits two essential properties: it is both strongly 
anisotropic and inhomogencous. Anisotropy comes from gravity, while inhomogeneity re- 
sults from the presence of top and bottom horizontal boundaries in convection cells. As 
a consequence, inertial range scalings of convcctivc turbulence, when they can ever be 
observed, depend strongly on the vertical coordinate. 

Accordingly, asymptotic theories of turbulence constructed under the assumptions of 
homogeneity and isotropy may be partially irrelevant to understand the observed prop- 
erties of turbulent convection. Neglecting intermittency effects, the cl assical picture re - 
garding scaling laws in this flow is that Bolgiano-Obukhov turbulence feolgianol(ll959l) : 
lObukhovl l|l959j) . hereafter B059 theory) should be present on correlation len gths r larger 
than the so-called Bolgiano length. A dimensional estimate of this length ( Chilla et all 
Il993|) is given by 

Nu 1/2 d 

LB = (RaPry/4 ' (L1) 
where d is the depth of the convective layer, Nu is the Nusselt number, Ra is the Rayleigh 
number and Pr is the Prandtl n umber. Instead, homogeneous and isotropic Kolmogorov 
turbulence l|Kolmogorovl l|l94l|) . hereafter K41 theory) should be observed for r < Lb- 
In the Bolgiano-Obukhov subrange, a dominant balance between buoyancy forcing (for 
an unstably stratified layer) and the third-order structure function occurs, which modi- 
fies the turbulent energy cascade substantially in comparison to K41: while longitudinal 
velocity increment scalings (Au r (r)) ~ r 1 / 3 are predicted asymptotically in the inertial 
range for K41, they should follow (Au r (r)) ~ r 3 / 5 in the B059 regime. As far as tempera- 
ture increments are concerned, (A6(r)) ~ r 1 / 3 results from K41 (passive scalar scalings), 
while (A9(r)) ~ r 1 / 5 is expected from dimensional arguments in B059 theory. These 
different regimes should in principle be detected through the measurements of structure 
functions or spectra when turbulence is sufficiently developed. However, clear evidences 
of K41 or B059 scalings are still lacking in both numerical and experimental convec- 
tion, even at very high Rayleigh numbers. Mixed scalings, which are not compatible with 
dimensional analysis, have sometimes been detected (e. g. K41 scalings for the veloc- 
ity and BQ59 scalings for temperature in the same wave number range, see for instance 
IVerzicco fc Camussi d2003l) ). BQ59 scalings have been re ported in the bulk of experimen- 



tal convection cells ( Benzi et (i^|l994| Ghing et all\20(Y$) , and in a numerical experiment 
ijCalzavarini. Toschi fc Trh)ic'cionell2002(h close to the horizontal walls of the convection 
cell. Direct scaling laws could not be observed for reduced structure functions (structure 
functions computed at fixed depth for correlation vectors r lying in horizontal planes 
only) in the latter study, so that the derivation of scaling exponents b v lCalzavarini et all 
l|2002l) relies on Extended Self Similarity (ESS, see lBenzi etall Jl99^). Re sults obtained 
for various anisotropic flows, such as the channel flow l|Arad et aZJ!l999<J) . demonstrate 
that it is sometimes possible to observe scaling laws by plotting structure functions di- 
rectly, even at modest Reynolds numbers, provided that the complete comp utation and 
SO (3) decomposition of structure functions are performed: in the work of lArad et all 
l)l999al) . the difference between reduced structure functions and the isotropic (angular 
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mean) component of the complete structure functions is striking. T his spherical harmon- 
ics decomposition of structure functions has also been computed bv lBiferale et all l|2003h 
for homogeneous Rayleigh-Benard convection (HRB) in a regime where Lb is compa- 
rable to the box size, thus preventing any possibility of B059 scalings. The anisotropic 
scaling exponents reported for HRB are anomalous, which means that they do not match 
dimensional analysis predictions either. The understanding of inertial range scalings in 
convective turbulence and their very existence even at high Raylcigh numbers therefore 
remains a puzzling problem that still deserves important efforts. 

This paper aims at presenting a complete description of turbulent statistics in Boussi- 
nesq convection at Ra = 10 6 and Pr = 1 in order to test explicitly the assumptions of 
both K41 and B059 theories in such a flow and to try to clarify some of the previously 
mentioned problems regarding the occurrence of inertial range scaling laws in convective 
turbulence. The study shows that inhomogeneity, anisotropy and the presence of buoy- 
ancy forcing at all scales all contribute to the absence of scaling laws, and points out 
that very high Raylcigh numbers should be reached in order to be able to observe the 
definite signature of spectral scalings in turbulent convection. The presentation of the 
results is as follows. The SO (3) decomposition of statistical averages is first used in 
to derive generalized Kolmogorov and Yaglom equations including inhomogeneous and 
anisotropic terms, which describe respectively the production, transport and dissipation 
of velocity and temperature fluctuations. A detailed numerical study of these equations 
is then proposed using Direct Numerical Simulation (DNS) data. The DNS and data pro- 
cessing algorithms are first described in The various contributions to the equations 
are then presented in ^]at the center of the convection cell and off the central plane, with 
particular emphasis laid on the anisotropic and inhomogeneous effects inferred from the 
analysis. The consequences and perspectives offered by these results are finally discussed 
in fj] 



2. Theoretical considerations 

Under the assumptions of homogeneity and isotropy, iKolmogorovl <)l94lh derived his 
famous equation for the third-order structure function of the turbulent velocity field 

((Aur) 3 ) = -| (e) r + 6^ ((Au r ) 2 ) , (2.1) 

which leads to the well-known K41 scaling law in the inertial range of turbulence. To ob- 
tain it, one also has to assume that the forcing mechanism of turbulence is important on 
large scales only. There are however numerous flows, such as Rayleigh-Benard convection, 
where some of these assumptions possibly break down. A detailed study of scale-by-scale 
energy budgets in turbulent convection is therefore required to test explicitly the ac- 
curacy of these assumptions. Indeed, the velocity budget will differ significantly from 
equation (|2.1[1 if one or several of these hypothesis are violated. The objective of this sec- 
tion is to provide a derivation of generalized forms of Kolmogorov and Yaglom equations 
including anisotropic and inhomogeneous effects, in order to be able to perform com- 
parisons with equation l|2.1|l and with the homogeneous and is otropic version of Yaglom 
equation, for a g i ven c o nvect iv e flow. A procedure in spire d by iMonin fc Yaeloml l|l975|) . 
iLindbord l|l996[) . Iflm] l|l997(l . lAntonia et d\ l|l997|) and iHnl ll2002>l. is adopted here . 
The method has already bee n applied to channel flow turbulence l|Danaih^_ei_ai] ]2j)0]J). 
heate d decaying turbulence (|Danaila et aZ.I fl999^ and shear turbulence l|Casciola et al\ 
l20f)3j) . which, as Rayleigh-Benard convection, are anisotropic (and sometimes inhomo- 
geneous) flows. Inhomogeneity is investigated in a similar fashion as in the paper by 
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iDanaila et oil l|200ll) . The proc edure is supplemen ted by a spherical harmonics decom- 
position of structure functions l|Arad et «/Jll999Sl) which serves to quantify anisotropic 
effects. It should be noted that a study of structure function scalings in turbulent con- 
vection is available jYak hotHl992T) and that it has been used by several authors (e. g. 
iBenzi. Toschi fc Tripicciond i|1998j) ') to predict scaling behaviour in numerical experi- 
ments. However, this work assumes isotropy and homogeneity, which, as will be shown 
later in the paper, are not appropriate for this flow. 

2.1. Derivation of the inhomogeneous Kolmogorov and Yaglom equations 

Summation over repeated indices is assumed throughout the derivation. A Boussinesq 
fluid of constant density p, kinematic viscosity v and thermal diffusivity k is used. Gravity 
is denoted by g = —q e z , where z is the vertical coordinate. The Boussinesq equations 
l)Chandrasekharlll96l|) for the velocity field u, pressure field p and temperature field T 
(separated into a mean part (T) and fluctuations 9) read 

d t Ui + u k d k Ui = ~-d t p + ag8 S iz + vdlui , (2.2) 
P 

d t (T) + d t e + u k d k e + u z d z (t) = K d 2 k e + K d 2 z (t) , (2.3) 

d k u k = (2.4) 



at position x and 



d t (I*) + d t 6' + u' k d' k 6' + u' z d' z (T) = K d' k 2 8' + k8' z 2 (T) , (2.6) 

d' k u' k = (2.7) 
at position x' (primes denote quantities evaluated at x'). As is usual in this kind of 
analysis, x and x' are independent variables so that djcU^ = d' k Ui = 0. Subtraction of 
equations l|2.2|l from (|2.5|) and of (|2.3|) from (|2.(i|l thus leads to 

<9 t AM 4 + Au k d k Au t + u k (d k + d' k ) Au. t = 

-- (di + d'i) Ap + ag A9 6 iz + v(d 2 k + d' k 2 ) A Ui , (2.8) 
P 

d t A (T) + t A6 + Au k d' k A9 + u k (d k + d' k ) AO + u' z d' z (T) - u z d z (T) = 

n(d 2 + d' 2 )Ae + K (d' 2 (T')^d 2 (T)) , 

(2.9) 

where Af = f — f for any variable /. Multiplying equations l|2.8|l and l|2.9() by 2Aui 
and 2A8 respectively, using incomprcssibility and averaging gives 



d t ((An,) 2 ) + (8' k {Au k {A Ul ) 2 )) + ((d k + d' k ) (u k (Au % ) 2 )) 

(A Ul ) 2 )-2((e) + (s 

(2.10) 



■- {(d t + di) (Au t A P )) + 2ag {A6Au z ) + v ({d 2 + d' 2 ) (Au t ) 2 ) - 2 ((e) + (e'» , 
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d t ((AO) 2 ) + (d' k (Au k (A9) 2 )) + ((d k + d' k ) (u k (A9) 2 )) 

+2 (u' z d' z (T 1 ) AO) - 2 (u z d z (T) AO) = 

K {(d 2 + d' 2 )(Ae) 2 )~2((N) + (N')) , (2.11) 



where (e) = v ((diUj) 2 ) and (N) = n((d t 6) 2 ). Following Iffill l|2002(K the analysis can 
be pushed further by using new variables X = (x + x')/2 and r = x' — x. Differential 
operators transform according to 



d 1 d 



dr k 2 dX k 

(2.12) 

& = — + i— 

k 8r k ^ 2 8X k - 

As the Rayleigh-Benard flow is homogeneous in the horizontal directions, the X depen- 
dence is restricted to Z, so that only dz derivatives have to be considered. It is important 
to outline that dz is not the same operator as d z , since r is to be kept constant when the 
former is applied, while z' must be kept constant in the latter case. Assuming statistical 
stationarity, equations H2.10fl - H2.llfl become 

A (Au k (A Ul ) 2 ) + \ A {Au z (Au t ) 2 ) + A ( Uz {Au t f) = 

-2 ((e) + (e')) - - A (Au z A P ) + 2ag (Au z A9) 
p oZ 

+2u£s ((A Ui ) 2 ) + \ A^ ((A Ui ) 2 ) , (2.13) 

A (Au k (A9f) + \ A <A«4A^) 2 ) + A ( Uz (Aef) 

+2 «^ (T') A0) - 2 (u 2 9, (T) A0) = 

-2 «JV) + (iV'}) + 2n A <(A0) 2 ) + f A_ ((A0)2) , 

(2.14) 

where all quantities now depend on (Z,r). 

2.2. Dealing with anisotropy 

To analyse the effects of anisotropy, a decomposition o f the r-dependence of statisti- 
cal averages on spherical harmonics must be performed l|Arad et aZ.lll999fc|) . For scalar 
averages, denoted here by (F), such as ((Au^) 2 ), this decomposition reads 

+oo +1 



(F)(Z,r)=J2 E {F) l m (Z,r)Yr(e^) 



=0 m=-£ 



where Y™(6, <p) denotes the spherical harmonic of sectorial order £ and azimuthal order 
to (here, 9 denotes the angle between the vertical direction z and r and should not 
be confused with a temperature perturbation). One also has to deal with the vectorial 
averages (U) = ((Au,) 2 Am) and (0) = ((A9) 2 Au) appearing in equations |2~l*3)l - l|2~l*4")l . 
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The spherical harmonics representations of these vectors are 



+00 +1 

(U)(Z,r) = J2 E (Un) e m (Z,r)RT(e^) 

+ (Us)t (Z, r)S?(0, <p) + (U T ) e m (Z, r)TT(9, <p) , 



(®)(z,r) = J2 E (eR)l(z,r)R?(e,v) 

t=0 m =-i 

+ (Gs)L (Z r)ST(e, + (©t)L (Z, r)TT((9, ¥>) , 

where = y £ m e r , SJ? 1 = VY e m , T v e l = V x flj 1 , and V is defined at r = 1. (C/ fl )L 
and (Us) l m are the components of a vector projected directly on a spherical harmonics 
vectorial basis. They are in that sense different from a quantity such as (Au z A6) , which 
results from a projection on e z followed by a projection on spherical harmonics. However, 
note that (Ur)^ could be noted ((Aui) 2 Au r ) m as well, since U™ is along e r . 

As the spherical harmonics form an orthogonal basis, equations for different I and m 
can be considered separately. Also, a dependence of statistical averages on Z and r only 
will be assumed from now on and the (Z, r) notation will be omitted subsequently. Using 
this formalism, the divergence operator acting for instance on (U) reads 



m ' 



while the laplacian of a scalar average is given by 



9 2 (F)V 18/j8 m i\ 1(1 + 1) , F s* 



' k 

Projecting equations (|2.13|) - (|2.14|l on Y™ for all I and m and applying the operator 

r 

1/+ 2 / y 2 • to the result, the following hierarchy of equations is readily obtained: 
Jo 



(ur)1 - / 



3 / r y 2 [(e) + ( £ ')]Ldy 
r Jo 

2 -^> j\*(Au z A9)lAy 



+ 2i ,-((A^) 2 > m <(A^) 2 ) m dy 



(2.15) 
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( e R/ m - 2 — / y(®s) m dy 



(i 



+ - (T) A< - «^ <T') A6)i) dy 

(2.16) 

In these equations, the inhomogeneous terms have been separated using the following 
notations: 

(NH Au y m = -^J\^(Au z (A Ul fy m dy 

and 

^Ol) = -^p 2 ^<^(^) 2 >ld 2/j 

The sums of these inhomogeneous terms in each equation will be used later on to char- 
acterize the net effect of inhomogeneity. They will be denoted by 

(NH) e m = {NH v ) l m + (NH p ) e m + (NH u y m + (NH Au f m , 
(NH°t = (NHl) l m + (NHl)l + (NH%J m . 



The usual Kolmogorov equation (|2.1() for homogeneous and isotropic turbulence can be 
straightforwardly recovered from equation H2.15(l . Using incomprcssibility, and isotropy, 
one finds that (U r ) = l/(3r 3 )d r (r 4 ((Au r ) 3 )) and d r ((Ait*) 2 ) = l/r 3 d r (r 4 9 r ((Au r ) 2 )) 
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(e. n. Iffiiil feooan . Furthermore, for homogeneous and isotropic turbulence (neglecting 
intermittency effects), (e) = (e') and (N) = (N 1 ) are constants and 

4 f V 2 i(e) + (e')t dy = ^ (e) rS m 6 m0 , (2.17) 
r Jo 6 

1 f y 2 [(iV) + <#')]£, dy = (iV) rS m S m0 . (2.18) 

In the preceding equations, the y/4n factor comes from the normalization with respect 
to Y °. Equation 1)2.1(1 follows directly from integration of equation If2.15f) by imposing 
£ = m = (isotropy ) and suppress i ng the (NH) terms (homogeneity). As noted by 
several authors (e. g. lAntonia et all lfl997tn . there is a striking analogy between the 
generalized Kolmogorov equation (|2.15|) which deals with the turbulent transport of 
kinetic energy (Ai^) 2 and equation (|2.16() . which involves the turbulent transport o f 
entropy fluctuations, which are proportional to {AO) 2 for a Boussinesq fluid l)L'vovll99l|) . 
This analogy is perfect when temperature is a passive scalar, while for non-neutrally 
stratified flows, extra terms representing the coupling between temperature and velocity 
fields appear. 

Finally, it is worth recalling that (Au z A9) m stands for the projection of the scalar cor- 
relator ((Am ■ e z )A9) on Y™. This notation may give the impression that equation l|2.15fl 
is closed in terms of correlation functions all belo nging to the s a me (£ , m) order. This 
is not the case because, as explained in detail bv lBiferale et al\ l)2003l) . gravity breaks 
the spherical symmetry of the equations of motion, leading to a coupling between the £ 
component of Kolmogorov equation and the (i — 1), £ and (£+ 1) harmonic components 
of the cross correlator between temperature and velocity. Therefore, projecting {Au z A0) 
on Y™ only provides a convenient way to analyse the net effect of gravity on the cor- 
responding £ component of the third-order structure function and does not mean that 
foliation of the Rayleigh-Benard equations with respect to £ (that is the question of the 
very existence of distinct inertial scaling exponents for different £ components of cor- 
relators) is satisfied. Foli ation remains an open que stion for both Rayleigh-Benard and 
Navier-Stokes equations rale fc Procacci 3 120051) . 

2.3. Symmetries of the Boussinesq equations 

Statistical averages in equations H2.15fl ~ H2.16fl obey some symmetry properties which will 
be listed here before a detailed analysis is carried out. First, there is an azimuthal sym- 
metry which is due to isotropy in horizontal planes, so that only the m — components 
need to be considered here. The other important symmetry is the top-down symmetry 
of the Boussinesq equations. For z = 1/2, the following rules can be derived: 

(a) (Ur) vanishes for odd £, 

(b) the even I components of statistical averages containing odd powers of u z , like for 
instance (u z (Aui) 2 Y , vanish, while their Z-derivatives are nonzero, 

(c) the Z-derivatives of the odd £ components of the same statistical averages vanish. 

Consequently, only even £ have to be considered at the center of the cell. Except for the 
azimuthal symmetry, the previous rules can not be applied to z ^ 1/2, where the complete 
£ spectrum has to be computed. An important consequence of the previous observations 
is that inhomogeneous effects should be present at the center of the convection cell, and 
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that the isotropic component should be affected by these effec ts. This also occurs for 
the channel flow, which obeys the same top-down symmetry i Danaila et a/jEoOlt) as 
Boussincsq convection. 



3. Numerics 

The various terms in equations (J2.15|) - (J2.16|) have to be computed in order to shed some 
light on the dominant contributions to the third-order structure functions (U) and (©) 
and to test quantitatively the assumptions of both B059 and K4f theories for this flow. 
To this end, one can choose to resort to either laboratory or numerical experiments. As 
will be discussed shortly, the main drawback of numerics is that only moderate Rayleigh 
numbers can be achieved in comparison to laboratory experiments. However, many inter- 
esting benefits can be gained from using numerical simulations to study the properties of 
equations iJ2.15J) - iJ2.16J) . The most important ones are that complete two-points averages 
can be computed easily without invoking the Taylor hypothesis and assuming isotropy, 
and that every quantity such as p can be easily obtained from the simulation at various 
locations. This motivated the choice of using numerics in the present study. A description 
of the DNS and numerical procedures used in order to perform this numerical analysis 
is now presented. 

3.1. Code setup 

To investigate equations H2.15JI - H2.16J) . the Boussinesq equations (|2.2J) - (|2.3J) for velocity 
and temperature are solved for R a = 10 6 , Pr = 1 and an aspect ratio A = 5 with a DNS 
code initiated bv lCalifanol 1 19961) . design ed for aniso tropic flows. It makes use of a sixth- 



order compact finite differences scheme l)LeleHl992|) in the vertical z direction and of a 
spectral scheme with 2 /3 dcaliasing in the horizontal ones . Time-stepping is implemented 
via a third-order, fully explicit Runge-Kutta scheme 1 Demuren. Wilson fc Caroenterl 



l200l|) . A projection method is used to ensure that the velocity field remains divergence- 
free. This method involves a Poisson solver to compute the pressure field. For the chosen 
aspect ratio, 256 points are used in each horizontal direction and 128 points in the ver- 
tical. A smaller mesh size is required in the vertical to resolve the thermal boundary 
layers accurately. Fixed temperature and stress-free boundary conditions are imposed on 
both plates, while periodic boundary conditions are used in the horizontal directions. 
Note that these boundary conditions offer a higher supercriticality (for a given Rayleigh 
number) than rigid boundaries, so that this simulation at Ra — 10 6 is in all likelihood 
located at the lower end (with respect to Ra) of the soft turbulence regime range of 
convection ijHeslot. Castaing fc Libchaberlll98?|l . 

Distances are normalized with respect to d: z = corresponds to the bottom plate and 
z = 1 to the top one. The simulation starts from small random temperature perturba- 
tions of the initial linear temperature profile and lasts for approximately fifty turnover 
times of the fluid, where one turnover time, defined here as twice the time it would take 
for a fluid element with velocity equal to the r.m.s. vertical velocity in the central plane 
to travel from one boundary to the other, corresponds to the time interval required before 
two snapshots can be considered independent. After a short linear growth phase, an ener- 
getically steady turbulent state i s reached in approximately one turnover t i me. It should 
be noted that plume clus tering l|Cattaneo et aZjEoOlt lHartlep et a/Jl20o3 IParodi et all 
120041 iRincon et aZ]l2005l) is observed in moderate and large aspect ratio simulations, 
which can be identified w ith the mean wind observ ed in small aspect ratio convection 
cells, as shown recently bv lvan Reeuwiik et all ()2005l) . The slow evolution of the integral 
scale towards large scales that results from this effect docs not affect the results pre- 
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sented here, which concern structures on scales smaller than the layer depth, for which 
an energetically steady state is reached quickly. In the statistically steady state, the Nus- 
selt number computed from the data of the simulation is Nu ~ 14. The Bolgiano scale, 
following equation (|1.1|) . is Lb — 0.11 (with respect to d), but as will be shown later, this 
value does not correspond to the effective Bolgiano scale at the center of the convection 
cell. The Kolmogorov dissipation scale, r\ = 0.016, is adequately resolved with the chosen 
mesh sizes. At the center of the convection cell, the Reynolds number based on the r.m.s. 
velocity and the layer depth is Re ~ 320 and the Taylor Reynolds number is Re\ ~ 30. 



3.2. Averaging procedure and statistical convergence 

Statistical averages are computed for given horizontal planes z by measuring each re- 
quired quantity once every turnover time and at e very four points in each horizontal 
direction. This procedure has notably been used bv lArad et al\ l)l999o|) for the channel 
flow and makes uses of homogeneity in both horizontal directions. At this stage, one 
may wonder why only a moderate l y turb ulent simulation with Ra = 10 6 is presented 
here, whereas IVerzicco fc Camussil ()2003h managed to simulate convcctive flows up to 
Ra = 2x 10 11 . It is worth emphasizing that the present work aims at computing various 
two-point statistical averages, and that such computations require a very long integration 
of the equations and a large number of independent realisations of the flow to be collected 
in order to get a good convergence, especially for structure functions of odd order and 
large I components. This very large data and computing time requirement is the main 
reason why all numerical studies that deal with the computation of statistical quantities 
such as structure functions have to make an important compromise between resolution 
and integration time a nd are therefore cur rently limited to moderate Rayleigh/Reynolds 
numbers. For instance. lArad et al\ l|l999af) used 1 60 independent snaps hots of their flow 
at Re\ — 70 to compute structure functions, while ICasciola et all ( 20031 ) simulated shear 
flow turbulen ce at Re\ = 45 during; 59 00 s hearing times in order to perform their statis- 
tical analysis. ICalzavarini et al\ l(2002|) and lBiferale et al\ (|2003h stored respectively 400 
and 270 independent configurations of aspect ratio one convective turbulence at Rayleigh 
numbers of order 10 7 . Such parameter va lues all lead to turbu lent regimes comparable 
to that of this simulation (in the case of lCasciola et al\ ()2003h . the largest value for r 
is 55^, to be compared with 30r/ here). In the present study, convergence for all statis- 
tical quantities at correlation lengths smaller than 0.5 starts to be satisfactory (up to 
a few pcrcents) after approximately ten turnover times, which is only possible because 
one time snapshot of a simulation with aspect ratio A = 5 can be considered ideally as 
25 independent snapshots of an aspect ratio one simulation for scales smaller than 0.5. 
Therefore, convergence is in fact only obtained by collecting more than 200 independent 
"aspect ratio one" snapshots, as in the previously mentio ned studies. In comparison, the 
high-resolution simulations bv lVerzicco fc Camussil (|2003h with aspect ratio A = 1/2 last 
for at most 220 free-fall turnover times, which correspond to only 16 turnover times as 
defined heref . This is far from sufficient to obtain convergence for all statistical quanti- 
ties of interest and, unfortunately, performing such simulations for much longer times is 
currently prohibitive from the computing time point of view. 



| Snapshots separated by one free-fall turnover time in their study can not be considered 
independent. The free-fall turnover time in their simulation is approximately 14 times smaller 
than the turnover time based on the r.m.s. vertical velocity, since the latter is approximately 
0.07 times the free-fall velocity (see their figure 14). 
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Figure 1. Horizontal maps of temperature fluctuations at z — 0.5 (left) and z — 0.9 (right). 

3.3. SO (3) decomposition 

The spherical harmonics decomposition is computed in the following way: statistical 
averages are computed for various correlation length vectors r making different angles 8 
with the vertical direction, located on a Gauss-Lcgcndre grid containing 15 points. This 
computation is performed for 20 different values of the azimuthal angle tp. A Fourier 
transform on tp and a Legcndre transform on 8 follow, which lead to the amplitudes of 
the various (£, m) components at a given r. The procedure is repeated at various distances 
r to obtain the final results. Interpolation between grid points is the only trick that one 
has to resort to in order to compute anisotropic structure functions. It does not affect the 
results significantly except for very small correlation lengths comparable to the mesh size. 
The dz and d\ derivatives in the inhomogencous terms are calculated by differentiating 
results obtained at different z with fourth-order finite difference schemes. 



4. Results 

This section provides an exhaustive presentation of results obtained by computing the 
various terms in equations i|2.15|l - ll2.1(jll from DNS data. A preliminary analysis is first 
carried out in order to outline inhomogencous effects in the Raylcigh-Bcnard system. 

4.1. Inhomogeneity in turbulent Rayleigh-Benard convection 

Figure n represents snapshots of temperature perturbations in the midlayer and close to 
the surface in the statistically stationary state of the simulation. Obviously, temperature 
maps at different altitudes look very different: while large scale temperature perturbations 
dominate the convective pattern inside the cell, smaller scale structures separated by 
dark lanes become prominent in the thermal boundary layer close to the surface. The 
presence of these dark lanes also leads to a c lear increase of sm all-scale power in the 
temperature spectra computed close the plates l|Rincon et q^jEoOsj) . This depth- variation 
of the properties of convective turbulence has been known for a long time. It is even 
more strikin g in non-Boussinesq c ompressible convection , where the top-down symmetry 
is lost (e. q. lToomre et al\ ljl990|) : ICattaneo et all l|l99 ll) ). Inhomogeneity in the system 
is further illustrated in figure [3 which shows horizontal velocity spectra computed at 
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Figure 2. Velocity spectra defined in equation 14. II for z — ( ), z = 0.13 ( ), 

2 = 0.25 (-■-■-■)> z = 0.38 ( — ■ ■ ■ — ). 2 = corresponds to the bottom of the convective cell. 
Due to the top-down symmetry of Boussinesq convection, only spectra in the lower half of the 
layer are presented. B059 and K41 slopes have been overplotted to illustrate the difficulty to 
identify a definite spectral slope. Wave numbers have been normalized with respect to 2k /(dA), 
so that k = 1 corresponds to the largest horizontal scale of the domain. 



various altitudes z. These spectra are defined according to 

E{k,z)= [ \u k {z)\ 2 kdn k , (4.1) 

where Q k stands for the wave vectors angles in the horizontal spectral plane, and the 
hat denotes horizontal Fourier transforms. Measuring definite slopes on these spectra 
proves all the more difficult than their shape depend quite strongly on the depth at 
which the computation is carried out. As will be demonstrated in the next paragraphs, 
there are various reasons for this absence of scaling laws. Note finally that the previously 
mentioned small-scale structures that appear on temperature maps close to the surface 
are essentially temperature perturbations, which explains why they do not show up in 
the corresponding velocity spectra at intermediate wave numbers. 



4.2. Isotropic component of equations \2.15)) - \2.1b)) 

The first step in the analysis of the simulation is the examination of the isotropic parts 
of equations (|2.15|l - l|2.16|l . which take on the following form: 

- WE + f y 2 ( A ^ A < <H/ + 2^ ((A Ui ) 2 ); + (NH)° = 

1 f >[(e) + ( £ ')]ody, (4.2) 
r Jo 
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Figure 3. Two different representations of the scale- by-scale budget equation 14. 21 . Left: lin- 

ear-linear plots of (U R )° ( ), -2ag/r 2 £ y 2 (Au z A9)l dy ( ), 2vd r ((Au,) 2 )" (- ■ - ■ 

-■), (NH}° ( ), -4/3 (e)r term ( ) and net budget ( ). Right: log-log plots 

of - (Ur)° ( ), 2ag/r 2 g y 2 <A U ,A0>° dy ( ), 2ud r ((Au,) 2 )° ( ), (NH)° (- ■ 

■ ■ — ). The sum of the l.h.s. terms of equation 14,21 ( ) is seen to be in very good balance 

with the r.h.s. 4/3 (e) r term ( ). 

4 [ r y 2 [(N) + (N')}° dy (4.3) 



o 

In the following, the right hand sides of equations H4.2[) - l|4.3[) wm be shown to be always 
close to 4\/47r/3 (e) r and 4\/47r/3 (N) r respectively, except on large scales r which cor- 
respond to penetration of the correlation length vector into the boundary layers, where 
dissipation (and so (e')) increases. For clarity reasons, these terms will thus be referred 
to as "4/3 (e) r" and "4/3 (N) r" terms. 

4.2.1. Central plane analysis 

Focus is first given to the central plane z = 1/2. An illustration of equation (|4.2|) at 
z = 1/2 is given in figure |3J The right hand side of equation l|4.2|l is almost linear with 
a slope 4\/47r/3 (e). This 4/3 (e) r term is the dominant contribution to the third-order 
structure function at moderate to large scales, but forcing by buoyancy is als o impo rtant 
in the same range, as seen on the left plot. As shown bv ICalzavarini et all l(2002t . the 
estimate of the Bolgiano length Hl.l(l is rather crude and the effective Bolgiano length 
(corresponding to the crossover between the buoyancy term and —4/3 (e) r) at the center 
of the convection cell which can be extrapolated from figure [3] is more likely closer to 1 
(r/rj ~ 60) than to 0.1 (r/rj = 6). The aspect ratio and set of boundary conditions have 
a direct influence on the prefactors of the various relations leading to equation i|l.l[l and 
thus affect the value of the effective Bolgiano length. Another simulation (not presented 
here) with aspect ratio 1 and all other parameters remaining constant notably reveals that 
Lb tends to get smaller in the central plane for convection cells with smaller aspect ratio. 
Th is effect may be rela ted to the presence of a mean wind (plume clustering): as shown 
bv lRincon et all 1 2005^ . the scale of the mean wind corresponds to the scale at which 



a maximum power is injected through buoyancy. In an aspect ratio 1 simulation, these 
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structures have to concentrate in a cube of aspect ratio 1, while they are allowed to spread 
over larger horizontal scales at larger aspect ratio. Buoyancy forcing may consequently 
be more important at smaller scales in small aspect ratio convection, thus reducing the 
effective Bolgiano length. Note finally that equation i|l.lfl does not take into account the 
fact that the Bolgiano length derived by Bolgiano and Obukhov depends on z-dependent 
dissipation rates and is therefore z-depe ndent as well. Using this exact expression, as 
done notably bv ICalzavarini et all l|2002() . does not lead to order of magnitude changes 
for Lb in the bulk of the convection cell in comparison to the estimation provided by 
equation 11. 1|) for the present simulation. In fact, both expressions rely on dimensional 
analysis and must therefore be understood as estimations up to order one prefactors. 

Using the effective value Lb — 1 for this simulation, one can estimate that Ls/r] ~ 60, 
so that there should in principle be some room for K41 scalings to occur for rj < r < Lb- 
However, the buoyancy term always remains of the same order of magnitude as the 
—4/3 (s) r term. Also, these two terms have opposite sign, which leads to a significant 
modification of the third-order structure function in comparison to isotropic turbulence 
forced on the largest scales only. The inhomogeneous term (NH)^ also gives a positive 
contribution that compensates the —4/3 (e) r term. This inhomogeneous term increases 
with r and becomes larger than the viscous term for r ~ 20?7. A detailed analysis of 
this term (figure 0| shows that the main contribution comes from (NH U X- The other 
inhomogeneous term s have a smaller influence on the global budget of equation (|4.2|) . 
iDanaila et al\ l|200lf) have shown that the (NH Uz ) term, which they have identified as a 
large-scale inhomogeneity, is also the dominant inhomogeneous contribution in channel 
flow turbulence. 

The most important consequence of having several important contributions to the 
third-order structure function is that it does not exhibit any scaling behaviour. To illus- 
trate this, a plot of (Ur)q is given in figure [3] together with a plot of the associated local 
logarithmic slope. It proves impossible to identify an inertial range on this plot, but note 
as a verification that the local slope of (Ur)q tends to 3 at very small scales, as expected 
from classical arguments. It is often claimed that the absence of scaling behaviour in 
various flows is due to the smallness of the Reynolds number, because in that case the 
viscous term in equation 1|2.1[) can not be neglected a most scales. This is only partly 
true for this convective flow, since the viscous term appears to be far smaller than the 
—4/3 (e) r term for r > 0.15 (r/r] > 10). Strictly speaking, the smallness of the Reynolds 
number is mainly responsible for the absence of a very large scale separation between 
Lb and rj here, which leads to the previously mentioned compensation effect. As will 
be discussed in jJSI definite scalings in the range 77 < r < Lb may only be observed for 
very large Raylcigh numbers. In any case, in a homogeneous flow forced on large scales 
only, for which no compensation by a forcing mechanism exist, scalings should in prin- 
ciple be observed, at least in a restricted subrange. There is however a supplementary 
effect which may prevent scaling laws from showing up even in this favourable situation, 
namely anisotropy. In that case, the occurrence of scalings may depend on the way data 
is analysed. This point will be further investigated later on in the paper. 

Comparing convective turbulence to other types of flows such as channel flows and 
shear flows shows that many features of turbulence are shared by all these flows, even 
though their forcing mechanisms are completely di fferent. For ins t ance, the right plot of 
figure |2| can easily be compared to the figure 5 of lDanaila et al\ l|200l|) . These authors 
studied the effect of inhomogeneity on third-order structure functions for the channel 
flow by assuming isotropy, which is however not completely satisfied in the flow. The 
main difference with the channel flow is that the (NH} term does not tend to balance 
—4/3 (e) r directly in turbulent convection at the largest scales that can be analysed us- 
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Figure 4. Decomposition of (NH)° ( ) into (NH V )° ( ), (NH p )l ( ), 

(NH Uz )l ( ), {NHau z )1 ( )• Tne lar S e scale inhomogeneity (NH Uz }° dominates 

the overall inhomogeneous contribution to equation l|4.2l) . 
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Figure 5. | (Ur)® \ as a function of r/n, normalized with respect to its maximum over r. Inset: 
logarithmic slope of the same function. A K41 scaling would lead to (Ur)q ~ r (up to a sign 
change). 

ing the SO (3) technique, since buoyancy forcing becomes increasingly important at such 
scales. Instead, the inhomogeneous term tends to compensate the net effect of —4/3 (e) r 
and of buoyancy (note that the situation is slightly different from the previously men- 
tioned channel flow study, since the SO (3) analysis is restricted to scales smaller than 
the integral scale of the flow, so that (Ur)q does not tend to ze ro yet at the largest sc ales 
available). Another interesting comparison with the results of lCasciola et all l)2003[l for 
homogeneous shear flow (their figure 1) and the left plot of figure |21 can be made. Their 
results are restricted to £ = 0, since their scale-by-scale budget results from averaging 
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Figure 6. Representations of the temperature scale- by-scale budget governed by equation l|4.3|l . 

- -), 2/r 2 /; y\{uA (T) A0)° - {u' z d' z (T') A9)° )dy 



Left: 
(--- 



linear-linear plots of (Qr) 
--), 2*d r ((A6) 2 )l (-■- 
and net budget ( 



(NH l 



-4/3 (N) r term 



Right: log-log plots of — (Qfl)o (" 



2/r 2 / r y 2 ({u z d z (T) A6)° - {u' z d' z <T'> A9)° )dy ( ), 2k8 t {{A6f)l ( ), (NH B )° Q (- 

■ ■ ■ — ). The sum of the l.h.s. terms of equation I4.3B ( ) is seen to be in very good balance 

with the r.h.s. 4/3 (N) r term ( ). 



structure functions over a r-sphere. Except for inhomogeneous effects, the scale-by-scale 
budgets of the two flows are very similar. Buoyancy appears to play almost the same role 
in Raylcigh-Bcnard convection as shear forcing in a shear flow. The same conclusions re- 
garding the determination of structure functions exponents can be drawn for both flows, 
since the forcing remains important at all scales in both cases. 

Equation 1)4. 3|1 for temperature fluctuations is also worth analysing in the central plane. 
Figure represents the various terms in this equation in a similar way as was done in 
figure |3 A plot of (Oij)o is presented in figure [7| together with the associated logarith- 
mic local slope. A clear plateau close to value 1 is observed for the logarithmic slope 
for 12?7 < r < 257/. This plateau corresponds to a dominant balance between (0_r)q and 

the —4/3 (N) r term, which is required in both B059 (V ~ ^r 3 /' 5 x (r 1 / 5 ) 2 ^ and K41 

~ ^V 1 / 3 x ( r1 ^ 3 ) 2 ^ theories. Figure El shows notably that the B059 hypothesis which 
neglects advection of the mean temperature profile term in comparison to the —4/3 (N) r 
term is very well satisfied in the simulation, because the core of the convection cell is 
almost isothermal at Ra = 10 6 , so that the d z (T) terms are negligible in this region. 
Inhomogeneity is seen to give a significant contribution to the third-order structure func- 
tion, mainly through the action of large-scale inhomogencitics (NH® r } (figure |HJ|, as 
in the velocity equation. However, this does not prevent scaling behaviour in a small 
"inertial" range in that case. The reason for this is that the dominant inhomogeneous 
term (NH® z )^, as can be seen on the same plot, is also proportional to r on large scales 
(right plot of figure El). Note that the observation of scalings for temperature related 
structure functions but not for velocity structure functions can not simply be attributed 
to the present value of the Prandtl number, because the potential scaling ranges of (Ur) 
and (Qr)q should be similar for Pr = 1. This analysis therefore confirms that moderate 
values of the Taylor Reynolds or Pcclct numbers (the equivalent of the Reynolds number 
for thermal diffusion) do not necessarily imply that an inertial range should not exist. 
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Figure 7. | (6_r)o I normalized to its maximum over r and its local logarithmic slope (inset), 

as a function of r/rj. 
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Figure 8. Decomposition of (NH B ) ( ) into (NH°) Q ( ), (NH* X ) Q ( ), 

(NH% U _^ ( ). As for the velocity equation, the large-scale inhomogeneity (NH® z ^ dom- 
inates the overall inhomogeneous contribution to equation l|4.3[l . 



4.2.2. Analysis off the central plane 

The main problem with the SO (3) decomposition of structure functions is that the 
analysis must be restricted to rather small scales r in regions close to the boundaries 
of the system (for wall- type boundary conditions). In such regions, the computation of 
reduced structure funct ions with r lying in plane s parallel to the boundaries remains the 
only available tool (see ICalzavarini et all (|2002|) for an application to Rayleigh-Benard 
convection). Remembering that the symmetry properties listed in ^disappear off the 
central plane, it is however interesting to push the SO (3) analysis as far as possible from 




Figure 10. Two representations of equation 14.31 for z = 1/4. The legends are the same as in 

figure El 



this region in order to observe trends of the z-dependence of inhomogeneous effects. A 
brief description of turbulent statistics at z = 1/4 is given to this end. At this altitude, the 
maximum r for which turbulent statistics can be computed is only I577, so that care must 
be taken in the interpretation of the results. For such a restricted range of correlation 
lengths, it is clear that at most a beginning of incrtial range can be expected and that 
probing inhomogeneous effects proves difficult, since these effects should be dominant on 
scales comparable to the integral scale of the flow only. 

As can be seen on figure and figure the behaviour of the various terms in equa- 
tions H4.2(l - (|4.3(l for r < lbrj and z = 1/4 is not very different from the case z = 1/2 
at first sight. There is however a slight difference concerning inhomogeneous effects in 
the velocity equation, which may be due to the closeness of the bottom wall: figure ITU 
reveals that normal pressure diffusion, which was small at small scales for z = 1/2, gives 
a significant contribution at the same scales for 2=1/4. This term is notably seen to 
continue to increase on the largest scales available, while (NH Uz ) seems to have reached 
a maximum for r close to the distance to the bottom boundary z. As far as the scaling 
of (C/r)q is concerned, figure fT51 shows a possible trend towards slope 1 for r/rj ~ 14, 
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Figure 12. | (Ur)q | normalized with respect to its maximum over r, as a function of r/n, at 
z — 1/4. Inset: logarithmic slope of the same function. 

but no real balance between the —4/3 (e) r term and the structure function is visible on 
figure [5J 

Meanwhile, figure ITTfl reveals that the logarithmic slope of (0_r)q, unlike for z = 1/2, 
does not seem to tend to 1 at the largest scales available for z = 1/4. Figurc^Jshows that 
inhomogcneous effects in the temperature equation are still mainly due to the (NH^ } 
term in this region. It therefore seems difficult to make definite conclusions regarding 
scaling laws in this region. It can only be asserted, following figure^ that the assumptions 
of homogeneity and of large-scale forcing are still clearly violated at z = 1/4. Explaining 
the logarithmic slope 1 for (Ur) q by simple K41 arguments therefore proves difficult in 
that case. This special slope may for instance be related to the penetration of r in the 
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Figure 13. | (©h)q | normalized with respect to its maximum over r and its local logarithmic 
slope (inset), as a function of r/rj, at z = 1/4. 




lower boundary layers where the flow dynamics are quite different from those in the bulk 
of the convection cell. 

4.3. Anisotropic sectors 

In the preceding paragraphs, two important assumptions of the K41 theory have been 
tested in the case of turbulent Rayleigh-Benard convection. The homogeneity hypothe- 
sis has been shown to be violated for r/rj > 20. The assumption of large-scale forcing 
has been shown to be very crude, since buoyancy forcing remains important at all non- 
dissipative scales of the simulation, even though the effective Bolgiano length is com- 
parable to the depth of the convective layer. It has thus been argued that the absence 
of scaling behaviour in turbulent convection at Ra = 10 6 was mostly related to these 
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effects, and that in the ideal situation of a homogeneous flow at Re\ = 30 forced on large 
scales only, inertial range scalings should be observed. There is however a possible supple- 
mentary cause for the apparent lack of iner tial range scalings i n turbulent flows, namely 
anisotropy. This has been demonstrated bv lArad et qj\ i|1999q|) for the channel flow and 
will be demonstrated in this paragraph for Rayleigh-Bcnard convection. An important 
remark is that in mildly turbulent numerical simulations, anisotropy may be caused by 
different effects. It may be due to physical effects such as spherical symmetry breaking 
by gravity or the presence of horizontal plates in the system, as outlined in the intro- 
duction, but may also be related to the fact that simulations are done in cartesian boxes 
th at break the rotational symmetry to o. This phenomenon has been notably pointed out 
bv lBiferale. Gualtieri fc Toschil who have shown that the projections of structure 

functions on anisotropic sectors were non- vanishing for a DNS of homogeneous isotropic 
turbulence performed in a cartesian domain. These various effects are difficult to disen- 
tangle. In order to quantify the exact contribution of spurious numerical anisotropy, it 
would be necessary to perform statistical analyses at much higher resolutions which, as 
mentioned in the section about numerics, is unfortunately not feasible currently. In every 
numerical study of anisotropic effects in turbulence, one should therefore keep in mind 
that spurious effects are present and consider the results with some caution. 

The analysis is restricted to the central plane z = 1/2. The largest sectorial order I 
considered in this paragraph is 6, for statistical convergence reasons. The scale-by-scale 
budgets for the anisotropic sectors read 

ITT \* , ITT \ l i 2a 9 r 2 /A AMii 

- (U R ) + (U s ) + — V (Au z A0) o dy 

4 f V [(c) + d» (4.4) 
r Jo 

and 

- (©h)o + (Qs)l + ^[v 2 < T > A < - «9' z (T>) A<) dy 

-I Ty 2 [{N) + (N')} e dy. 
r Jo 

(4.5) 

These budgets involve extra terms related to anisotropy in comparison to equations Il4.2|l - 
H4.3(l . The most important ones are those containing the transverse third-order structure 
functions (Us)q and (Os)q- To illustrate this, the I = 2 component of equations l|2.15|l - 
H2.16fl can be investigated in more detail Cfigure H3)l . The (positive) transverse third-order 
structure function terms are shown to compensate significantly the (negative) longitudi- 
nal third-order structure functions in each equation. There is also a small anisotropic con- 
tribution of the [(e) + (e')]o an< ^ + (^')]o terms at the largest r, which is due to in- 
creased dissipation in the boundary layers. Finally, as was the case for the isotropic equa- 
tions, the inhomogeneous terms contribute to the scale-by-scalc budgets in the anisotropic 
sectors. 

Figure El displays the spherical harmonics spectrum of (Ur) 1 for even I up to I = 6 
(remember that odd £ contributions are zero when considering averages in the central 
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Figure 15. Left: equation (EHSt for £ = 2 with ([/a)* ( ), -£{£ + l)/r 2 y (U s )' dy 

(+++++), — 2/r 2 j/ 2 [(e) + (e')] dy ( ), buoyancy term ( ), viscous term (-■-■- 

■), (NH)l ( ) and net budget ( ). Right: equation fTTGl for £ = 2 with (Q R ) e (— 

- -), -£(£ + l)/r 2 / r y {@ s }* dy (+++++), -2/r 2 gy 2 \(N) + {N% dy ( ), mean tern- 
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Figure 16. Spherical harmonics spectrum of (Ur) at the center of the convection cell for 

r/77 = 7.9 ( )j r/rj = 15.3 ( ) and r/77 = 22.6 ( ). Spectra are normalized 

with respect to (Ur)q (r) for each r. The 1 — 2 and £ = 4 terms are shown to have a comparable 
magnitude to the £ — one at moderate to large correlation lengths. 



plane). It is once again emphasized that this spectrum contains contributions of both 
physical anisotropics due to gravity and spurious numerical anisotropies due to the use 
of a cartesian computational domain. The I 7^ components appear to be significant in 
comparison to their isotropic counterpart at moderate to large correlation lengths. As the 
complete third-order structure functions result from a linear combination (with weights 
given by the spherical harmonics functions) of these various significant contributions, 
which may follow distinct scaling laws if foliation occurs, it is expected that genuine 
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incrtial range scalings predicted by isotropic theories should be difficult to observe for 
reduced structure functions, which correspond to the complete structure functions taken 
at a polar angle 9 = tt/2 only and therefore involve the contribution of all even spherical 
harmonics components of these complete structure functions evaluated at this angle. 
This is confirmed by figure El an d figure El figure El shows a reconstruction of the 
reduced structure functions (U r ) (9 = tt/2) and (0 r ) (9 = tt/2) averaged over using the 
£ = 0,2,4,6 components of the complete (U r ) = ((A-u i ) 2 Aw r ) and (6 r ) = ((A#) 2 Au r ) 
for the present simulation. Figure 1181 is a similar reconstruction of the local logarithmic 
slopes of the same objects. The shape and local logarithmic slopes of the reduced structure 
functions differ substantially from those of the £ = components, which shows that 
anisotropic components (especially £ = 2, as can be seen on both plots of figure I17J) 
have a marked influence on the local scaling exponents of reduced structure functions. In 
the case of (Or), one finds for instance that the local exponent of the reduced structure 
function for r/r\ = 20 is approximately 0.8, while the £ — component exponent is 1 
at the same correlation length, as shown earlier in the paper and on the right plot of 
figure IT51 A similar analysis stands for (U r ) (left plot of figure ITB|l . Such differences can 
not be ignored, since they are larger than typical intermittency corrections to inertial 
range scalings and than differences between concurrent theories of turbulence. Inertial 
range scaling exponents obtained via an analysis of reduced structure functions should 
therefore be regarded with caution, at least when the flow is fully three-dimensional. For 
the same reasons, computations of structure functions in laboratory experiments, which 
rely on single points measurements and on the Taylor hypothesis along a specific direction 
(that of the mean wind in Rayleigh-Benard convection) will not necessarily lead to correct 
measurements of scaling exponents predicted by isotropic theories or to reliable values of 
intermittency corrections, because such exponents only apply to the isotropic component 
of the structure functions. Note however that the previous conclusions do not mean that 
measuring structure functions along a specific direction does not make physical sense. 
For instance, using reduced structure functions may prove useful to analyse scalings in 
two-dimensional flows, layered flows or in strongly anisotropic flows such as boundary 
layers, where scalings in the directions perpendicular to the shear can be observed. Also, 
in regions where the SO (3) decomposition can not be performed, as is the case close to 
walls, the only currently available methods of analysis are the computation of reduced 
structure functions and single points measurements coupled to the Taylor hypothesis. 



The present analysis finally proves useful to test theoretical arguments on anisotropic 
turbulence and the phenomenon of foliation: a ccording to modern theories of anisotropy 
l)Arad et oi.lll999ll iBiferale fc Procacciall2005h , each £ component of structure functions 
may follow a distinct inertial range scaling law with a scaling exponent increasing with £ . 
Although it proves difficult to identify scaling ranges (and therefore scaling exponents) 
on the £ ^ sectors of the third-order structure functions in this simulation, the increase 
of scaling exponents with £ seems to be qualitatively well verified, at least it is consistent 
with the fact that the £ = 2 contribution is more important for large r than for small r 
in comparison to the £ = component. This may be viewed as an argument in favour of 
a small-scale return to isotropy. Note however that definite conclusions can not be drawn 
here because the £ = 2 component still remains significant for small r (figure IT?}) . This 
may be a finite resolution effect related to the use of cartesian grids and boxes that do 
not have rotational symmetry 
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Figure 17. Reconstruction of the reduced structure functions ([/, ) (6 = tt/2) (left plot, + 
symbols) and (Or) (0 = tt/2) (right plot, + symbols) at the center of the convection cell using 

anisotropic components up to £ = ( ), £ — 2 ( ), £ = 4 (-■-■-■), f = 6 ( — ■ ■ ■ — ). 

Note the significant contribution of the £ = 2 components to the reduced structure functions. 




Figure 18. Reconstruction of the local logarithmic slope of reduced structure functions 

(U r ) (0 = 7r/2) (left plot, ) and (9 r ) (0 = tt/2) (right plot, ) at the center of the 

convection cell, using anisotropic components up to £ — (+), 1 = 2 (x), 1 = 4 (o), I = 6 (A). 



5. Summary and discussion 

New results on inertial range scaling laws, scalc-by-scale budgets and correlation func- 
tions in turbulent Rayleigh-Benard convection at moderate Raylcigh and Reynolds num- 
bers have been reported in this paper. A derivation of generalized Kolmogorov and Ya- 
glom equations 12.15fl - 12.16f) . which notably makes use of the SO (3) decomposition of 
statistical averages, has first been presented. The formalism provides a convenient way of 
disentangling inhomogeneous and anisotropic effects in convective turbulence, but it can 
also be used to study other anisotropic and/or inhomogeneous numerically simulated or 
experimental flows. For this purpose, the buoyancy forcing term in the velocity equation 
should be replaced by any appropriate anisotropic forcing mechanism. 

The analysis of scale-by-scale budgets in a convection DNS at Ra = 10 6 has led to the 
following conclusions. First of all, a compensation effect between the classical —4/3 (e) r 
term of the generalized Kolmogorov equation on one hand, inhomogeneous terms, viscous 
terms and buoyancy forcing on the other hand, has been demonstrated in the isotropic 
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sector of the velocity equation. Among these three terms, buoyancy is the most impor- 
tant one. This compensation effect prevents any clear scaling law from showing up at 
Re\ = 30 in the case of the third-order velocity structure function. Also, equation i|l.lfl 
has been shown to underestimate significantly the effective value of the Bolgiano length, 
at least for convective flows in moderate aspect ratio containers. The analysis of the 
mixed temperature structure function (0_r)q has revealed a scaling behaviour compat- 
ible with both K4f and B059 theories in the central plane, in spite of the presence of 
inhomogeneous effects on large scales. Meanwhile, neglecting the term related to advec- 
tion of the mean temperature profile, which is required in the B059 theory, has been 
shown to be justified. Scale-by-scale budgets computed off the central plane have led to 
similar conclusions. An interesting inhomogeneous effect at z = 1/4 is the significant 
contribution of pressure diffusion normal to the bottom wall. 

A spherical harmonics spectrum of structure functions has finally revealed that the low 
£ degree contributions to the third-order structure functions are not small in comparison 
to their isotropic counterparts at moderate to large scales. It has particularly been shown 
that using reduced structure functions at fixed polar angle in the bulk of a convection cell 
in order to reveal scaling behaviour predicted by isotropic theories is misleading, since 
these structure functions involve linear combinations of various I components which 
scale differently with respect to r. This argument also indicates that using single points 
measurements in the bulk of three-dimensional flows together with the Taylor hypothesis 
in the particular direction of a mean flow to test the predictions of asymptotic dimensional 
isotropic theories of turbulence or to calculate intermittency corrections to these theories 
may lead to significant biases. It has however been outlined that using reduced structure 
functions may still be appropriate when the direction along which these functions are 
computed has a special importance, as is the case for instance in layered flows or boundary 
layers. Also, it has been pointed out that spurious anisotropics related to the use of 
cartesian domains in simulations of mildly turbulent flows, which are difficult to quantify, 
certainly interfere with the genuine anisotropic effects due to gravity. Thus, the respective 
amplitudes of the various £ components of structure functions at moderate Reynolds 
number has to be considered with some caution. Independ ently of the problem of the 
origin of anisotropies, these results confirm the analysis of lArad et al\ l)l999oj) for the 
channel flow: disentangling anisotropic effects can not be avoided if one wants to study 
the scaling behaviour of anisotropic systems correctly. Besides, even though the study of 
the r-dependence of different i components in equations H2.15JI - H2.16J) seems to indicate 
qualitatively small-scale return to isotropy in the case of Raylcigh-Bcnard convection, 
it can not be ruled out that anisotropic effects persist at small scales, even at very 
high Reynolds numbers. Small-scale anisotropies (called "ramp-cliff struc tures" ) have for 
instan c e been reported in studies o f turbule nt mixing of pas sive scalars ( Mestaver et all 
l|l976h : lAntonia & Van Attal <|l978h . see also IWarhaftl Ij200(t) ). 

An important issue is to which extent the present results can be used to infer scaling 
behaviour at very high Rayleigh and Reynolds numbers. The still important contribution 
of buoyancy forcing to the scale- by-scale budget at scales smaller than Lb naturally 
raises the question of how many scale decades are required in a turbulent convective 
flow in order to be able to identi fy definite scalings. In the s oft turbulence regime, where 
Nu ~ Ra? with 1/4 < 7 < 1/3 l(Grossmann & Lohsell27)0tl . the Bolgiano length should 
not depend significantly o n Ra. For instance, follow ing equation (Jl.lf) . Lb ~ i?a -3 / 28 
is obtained for 7 = 2/7 IJProcaccia fc Zeitakl 119891) . In the hard regime where it has 
recently been shown that Nu ~ (RaPr) 1 / 2 ~ ijCalzavarini et all\200^i , the Bolgiano length 
should remain constant (with respect to both Ra and Pr in that case, while the Pr 
dependence in the hard turbulence regime depends on whether Pr > 1 or Pr < 1), 
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thus leading to the same result at all Rayleigh numbers. For the present aspect ratio 
and Prandtl number, Lb should therefore remain close to 1 at all Rayleigh numbers. 
On the contrary, the ratio between the dissipation scale and the effective Bolgiano scale 
is expected to increase with increasing Rayleigh number, w hich should in principle hel p 
to observe K41 scalings if rj <C r <C Lb can be achieved l|Grossmann fc L~ vl 1199.1) . 
One should however derive more precise conditions of applicability of K41 in that case. 
Assuming K41 to be valid for r < Lb, the buoyancy term, which equals the —4/3 (e) r 
term at Lb, would scale like r 5 / 3 in that range. This approximation looks rather crude 
for r close to Lb, but as the buoyancy term should scale like r 9 ' 5 above Lb according 
to B059, the actual local scaling exponent close to Lb may not be very different from 
this 5/3 value. The important point here is that this exponent be larger than 1 (which is 
the exponent of the —4/3 (e) r term). According to these scalings, non-dissipative scales 
as small as 3 x 10~ 2 Ls should be available in order for the buoyancy term to become 
ten times smaller than the —4/3 (e) r term, which roughly corresponds to the conditions 
of the present simulation, in which not definite K41 plateau can be observed (figure |3J). 
In order to find a range of scales in which the buoyancy term would become at least one 
hundred times smaller than the —4/3 (e) r term, one should then look for a regime where 
77 < 1CP 3 £ r. The required Rayleigh n umber can be evaluated by using 77 ~ i?a~ 9 / 28 , 
proposed by Grossmann fc Lohsel l)l993|) . This relation can be calibrated using the results 
of the present simulation at Ra = 10 6 , for which 77 = 0.016. One finally finds that Ra has 
to exceed several times 10 9 to obtain a flow with 77 < 10~ 3 Ls. Even in this situation, 
definite scalings are expected only in the subrange IO77 < r < O.ILb, that is over one 
decade. The main problem is that two scale separations r -C Lb and r 3> 77 must 
simultaneously be satisfied, which can only occur at very high Rayleigh number. This is 
an important argument to understand why inertial range scalings have been difficult to 
identify in experimental or numerical turbulent convection studies. 

For the preceding reasons, it is probabl e that various results on scaling exponents ob- 
tained in the range 10 6 < Ra< 10 9 (e. a. ICalzavarini et all l|2002l) . who have computed 
local scaling exponents from reduced structure functions) are biased simultaneously by 
anisotropic effects and significant buoyancy forcing at r < Lb- The lack of complete K41 
or BQ59 scalings in the sim ulations at very high Rayleigh number (up to Ra = 2 x 10 11 ) of 
IVerzicco fc Camussil l)2003|) may also be a hint that inertial range scalings in anisotropic 
flows such as Rayleigh-Benard convection can not be understood easily using isotropic 
theories. If inhomogencity and isotropy were to play an important role in the scaling 
behaviour of very high Rayleigh number convective turbulence, the present results could 
serve as guidelines to study them in detail. More generally, the SO (3) decomposition of 
structure functions appears to be a very powerful (and necessary) tool to study various 
anisotropic flows. It should be particularly helpful to use it when possible in research 
domains such as astrophysics or atmospheric sciences, where anisotropic turbulence is 
ubiquitous. Such an investigation for strongly stratified non-Boussinesq turbulent con- 
vection, which is relevant to the Sun and other stars, is currently underway. 



The author acknowledges several fruitful discussions with A. A. Schekochihin and 
thanks F. Ansclmet, F. Lignicrcs and M. Rieutord for many insightful comments and 
criticisms on early versions of the manuscript. Numerical simulations have been per- 
formed on the IBM SP4 supercomputer of Institut du Developpement et des Ressources 
en Informatique Scicntifiquc (IDRIS, Orsay, France), which is gratefully acknowledged. 



Anisotropy, inhomogeneity and inertial range scalings in turbulent convection 27 

REFERENCES 

Antonia, R. & Van Atta, C. 1978 Structure functions of temperature fluctuations in turbulent 
shear flows. J. Fluid Mech. 84, 561-580. 

Antonia, R. A., Ould-Rouis, M., Anselmet, F. & Zhu, Y. 1997 Analogy between predic- 
tions of Kolmogorov and Yaglom. J. Fluid Mech. 332, 395-409. 

Arad, I., Biferale, L., Mazzitelli, I. & PROCACCIA, I. 1999a Disentangling scaling prop- 
erties in anisotropic and inhomogeneous turbulence. Phys. Rev. Lett. 82, 5040-5043. 

Arad, I., L'vov, V. & PROCACCIA, I. 19996 Correlation functions in isotropic and anisotropic 
turbulence: the role of the symmetry group. Phys. Rev. E 59, 6753-6765. 

Benzi, R., Ciliberto, S., Tripiccione, R., Baudet, C, Massaioli, F. & Succi, S. 1993 
Extended self-similarity in turbulent flows. Phys. Rev. E 48, 29. 

Benzi, R., Toschi, F. & Tripiccione, R. 1998 On the heat transfer in Rayleigh-Benard 
systems. J. Stat. Phys. 93. 

Benzi, R., Tripiccione, R., Massaioli, F., Succi, S. & Ciliberto, S. 1994 On the scaling of 
the velocity and temperature structure functions in Rayleigh-Benard convection. Europhys. 
Lett. 25 (5), 341-346. 

Biferale, L., Calzavarini, E., Toschi, F. & Tripiccione, R. 2003 Universality of 
anisotropic fluctuations from numerical simulations of turbulent flows. Europhys. Lett. 64, 
461-467. 

Biferale, L., Gualtieri, P. & Toschi, F. 2000 Statistics of pressure and of pressure- velocity 

correlations in isotropic turbulence. Phys. Fluids 12 (7), 1836. 
Biferale, L. & Procaccia, I. 2005 Anisotropy in turbulent flows and in turbulent transport. 

Phys. Reports 414 (2), 43. 
Bolgiano, R. 1959 Structure of turbulence in stratified media. J. Geophys. Res. 64, 2226. 
Califano, F. 1996 A numerical algorithm for geophysical and astrophysical inhomogeneous 

fluid flows. Comp. Phys. Comm. 99, 29. 
Calzavarini, E., Lohse, D., Toschi, F. & Tripiccione, R. 2005 Rayleigh and Prandtl 

number scaling in the bulk of Rayleigh-Benard turbulence. Phys. Fluids 17 (5), 055107. 
Calzavarini, E., Toschi, F. & Tripiccione, R. 2002 Evidences of Bolgiano scaling in 3D 

Rayleigh-Benard convection. Phys. Rev. E 66, 016304. 
Casciola, C. M., Gualtieri, P., Benzi, R. & Piva, R. 2003 Scale-by-scale budget and 

similarity laws for shear turbulence. J. Fluid Mech. 476, 105-114. 
Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A. & Hurlburt, N. E. 1991 

Turbulent compressible convection. Astrophys. J. 370, 282-294. 
Cattaneo, F., Lenz, D. & Weiss, N. 2001 On the origin of the solar mesogranulation. Astro- 
phys. J. 563, L91-L94. 
Chandrasekhar, S. 1961 Hydrodynamic and hydromagnetic stability. Dover. 
Chilla, F., Ciliberto, S., Innocenti, C. & Pampaloni, E. 1993 Boundary layer and scaling 

properties in turbulent thermal convection. Nuovo Cimento 15D (9), 1229. 
Ching, E. S. C, Chui, K. W., Shang, X.-D., Qiu, X.-L., Tong, P. & Xia, K.-Q. 2004 

Velocity and temperature cross-scaling in turbulent thermal convection. J. Turbulence 5. 
Danaila, L., Anselmet, F., Zhou, T. & Antonia, R. A. 1999 A generalization of Yaglom 's 

equation which accounts for the large-scale forcing in heated decaying turbulence. 391 391, 

359-372. 

Danaila, L., Anselmet, F., Zhou, T. & Antonia, R. A. 2001 Turbulent energy scale budget 

equations in a fully developed channel flow. J. Fluid Mech. 430, 87-109. 
Demuren, A. O., Wilson, R. V. & Carpenter, M. 2001 Higher-order compact schemes for 

numerical simulation of incompressible flows, Part I: theoretical development. Numerical 

Heat Transfer 39, 207-230. 
GROSSMANN, S. & Lohse, D. 1993 Characteristic scales in Rayleigh-Benard turbulence. Phys. 

Lett. A 173, 58. 

GROSSMANN, S. & Lohse, D. 2000 Scaling in thermal convection: a unifying theory. J. Fluid 
Mech. 407, 27. 

GROSSMANN, S. & L'vov, V. 1993 Crossover of spectral scaling in thermal turbulence. Phys. 
Rev. E 47, 4161. 



28 



F. Rincon 



Hartlep, T., Tilgner, A. & Busse, F. H. 2003 Large scale structures in Rayleigh-Benard 
convection at high Rayleigh numbers. Phys. Rev. Lett. 91 (6), 064501. 

Heslot, F., Castaing, B. & Libchaber, A. 1987 Transitions to turbulence in helium gas. 
Phys. Rev. A 36, 5870-5873. 

Hill, R. J. 1997 Applicability of Kolmogorov's and Monin's equations of turbulence. J. Fluid 
Mech. 353, 67-81. 

Hill, R. J. 2002 The approach of turbulence to the locally homogeneous asymptote as studied 

using exact structure-function equations. arXiv.physics/0206034 ■ 
KOLMOGOROV, A. N. 1941 Dissipation of energy in locally isotropic turbulence. Dokl. Akad. 

Nauk SSSR 32, 16. 

Lele, S. K. 1992 Compact finite difference schemes with spectral-like resolution. J. Comp. 
Phys. 103(1), 16-42. 

LlNDBORG, E. 1996 A note on Kolmogorov's third-order structure- function law, the local 
isotropy hypothesis and the pressure-velocity correlation. J. Fluid Mech. 326, 343-356. 

L'vov, V. S. 1991 Spectra of velocity and temperature fluctuations with constant entropy flux 
of fully developed free-convective turbulence. Phys. Rev. Lett. 67, 687-690. 

Mestayer, P., Gibson, C. H., Coantic, M. & Patel, A. 1976 Local anisotropy in heated 
and cooled turbulent boundary layers. Phys. Fluids 19, 1279-1287. 

Monin, A. S. & Yaglom, A. M. 1975 Statistical fluid mechanics, , vol. 2. MIT Press. 

Obukhov, A. M. 1959 Dokl. Akad. Nauk. SSR 125, 1246. 

Parodi, A., von Hardenberg, J., Passoni, G., Provenzale, A. & Spiegel, E. A. 2004 
Clustering of Plumes in Turbulent Convection. Phys. Rev. Lett. 92 (19), 194503. 

PROCACCIA, I. & Zeitak, R. 1989 Scaling exponents in nonisotropic convective turbulence. 
Phys. Rev. Lett. 62, 2128. 

van Reeuwijk, M., Jonker, H. J. J. & Hanjalic, K. 2005 Identification of the wind in 
Rayleigh-Benard convection. Phys. Fluids 17, 051704. 

Rincon, F., Lignieres, F. & Rieutord, M. 2005 Mesoscale flows in large aspect ratio simu- 
lations of turbulent compressible convection. Astron. & Astrophys. 430, L57-L60. 

Toomre, J., Brummell, N., Cattaneo, F. & Hurlburt, N. E. 1990 Three-dimensional 
compressible convection at low Prandtl numbers. Comp. Phys. Comm. 59, 105-117. 

Verzicco, R. & Camussi, R. 2003 Numerical experiments on strongly turbulent thermal con- 
vection in a slender cylindrical cell. J. Fluid Mech. 477, 19-49. 

Warhaft, Z. 2000 Passive scalars in turbulent flows. Ann. Rev. Fluid Mech. 32, 203-240. 

Yakhot, V. 1992 4/5 Kolmogorov law for statistically stationary turbulence: application to 
high-Rayleigh-number Benard convection. Phys. Rev. Lett. 69, 769. 



